#!/usr/bin/perl
use strict;
use warnings;
#DESC: Determines best alignment between two sequences using the Needleman-Wunsch
#     algorithm given a file with two sequences and a file containing a scoring
#     matrix and gap score.
#USAGE: perl seqAlignment.pl [matrix filename] [sequence filename]

my @align = needlemanWunsch("adamwest","admwst");

for (my $i=0; $i<=length($align[0]); $i+=80){
  print substr($align[0],$i,80)."\n";
  print substr($align[1],$i,80)."\n";
  print "\n";
}

# the needleman-wunsch pairwise alignment algorithm, courtesy of the be131 homework solutions
sub needlemanWunsch {
  my (@seqs) = @_;
  my @s = ();
  my @tb = ();
  my $tempMax = 0;
  my $pair = '';
  my $nextI = 0;
  my $nextJ = 0;
  my $nextPos = '';
  my $seq1 = '';
  my $seq2 = '';
  my $flagEnd =0;

  for(my $i=0;$i<2;$i++) {
    $seqs[$i] = uc $seqs[$i];
  }

  my $gapPenalty = -4;
  my @matrixChars = ("A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V", "B", "Z", "X");
  my %preMatrix = 
    ( "A" => [4, -1, -2, -2, 0, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -3, -2, 0, -2, -1, 0],
      "R" => [-1, 5, 0, -2, -3, 1, 0, -2, 0, -3, -2, 2, -1, -3, -2, -1, -1, -3, -2, -3, -1, 0, -1],
      "N" => [-2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3, 0, -2, -3, -2, 1, 0, -4, -2, -3, 3, 0, -1],
      "D" => [-2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4, -3, -3, 4, 1, -1],
      "C" => [0, -3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2],
      "Q" => [-1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2, -1, -2, 0, 3, -1],
      "E" => [-1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3, 1, -2, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1],
      "G" => [0, -2, 0, -1, -3, -2, -2, 6, -2, -4, -4, -2, -3, -3, -2, 0, -2, -2, -3, -3, -1, -2, -1],
      "H" => [-2, 0, 1, -1, -3, 0, 0, -2, 8, -3, -3, -1, -2, -1, -2, -1, -2, -2, 2, -3, 0, 0, -1],
      "I" => [-1, -3, -3, -3, -1, -3, -3, -4, -3, 4, 2, -3, 1, 0, -3, -2, -1, -3, -1, 3, -3, -3, -1],
      "L" => [-1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, -2, 2, 0, -3, -2, -1, -2, -1, 1, -4, -3, -1],
      "K" => [-1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, -1, -3, -1, 0, -1, -3, -2, -2, 0, 1, -1],
      "M" => [-1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, 0, -2, -1, -1, -1, -1, 1, -3, -1, -1],
      "F" => [-2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, -4, -2, -2, 1, 3, -1, -3, -3, -1],
      "P" => [-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, -1, -1, -4, -3, -2, -2, -1, -2],
      "S" => [1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, 1, -3, -2, -2, 0, 0, 0],
      "T" => [0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, -2, -2, 0, -1, -1, 0],
      "W" => [-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, 2, -3, -4, -3, -2],
      "Y" => [-2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7, -1, -3, -2, -1],
      "V" => [0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4, -3, -2, -1],
      "B" => [-2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4, 1, -1],
      "Z" => [-1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1],
      "X" => [0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1]);

  my %simMatrix;
  for my $amino (keys %preMatrix) {
    for(my $i=0;$i<@matrixChars;$i++) {
      $simMatrix{$amino . $matrixChars[$i]} = $preMatrix{$amino}[$i];
    }
  }

#THIS SECTION PERFORMS THE NEEDLEMAN-WUNSCH ALGORITHM

#initialize the array with gap penalties
  $s[0][0] = 0;
  $tb[0][0] = 'END';
  for (my $i=1; $i<=length($seqs[1]); $i++){
    $s[$i][0] = $i * $gapPenalty;
    $tb[$i][0] = 'U';
  }
  for (my $j=1; $j<=length($seqs[0]); $j++){
    $s[0][$j] = $j * $gapPenalty;
    $tb[0][$j] = 'L';
  }

#perform alignment algorithm, store scores in array @s and traceback
#info in array @tb
  for (my $i=1; $i<=length($seqs[1]); $i++){
    for (my $j=1; $j<=length($seqs[0]); $j++){
      $pair = substr($seqs[1],$i-1,1) . substr($seqs[0],$j-1,1);
      $s[$i][$j] = $s[$i-1][$j-1] + $simMatrix{$pair};
      $tb[$i][$j] = 'D';
      $tempMax = $s[$i-1][$j] + $gapPenalty;
      if ($tempMax > $s[$i][$j]) {
        $s[$i][$j] = $tempMax;
        $tb[$i][$j] = 'U';
      }
      $tempMax = $s[$i][$j-1] + $gapPenalty;
      if ($tempMax > $s[$i][$j]) {
        $s[$i][$j] = $tempMax;
        $tb[$i][$j] = 'L';
      }
    }
  }

#perform traceback and print alignment
  $nextI = length($seqs[1]);
  $nextJ = length($seqs[0]);
  $nextPos = $tb[$nextI][$nextJ];
  while ($flagEnd == 0) {
    if ($nextPos eq 'D') {
      $seq1 = substr($seqs[0],$nextJ-1,1) . $seq1;
      $seq2 = substr($seqs[1],$nextI-1,1) . $seq2;
      $nextI--;
      $nextJ--;
    }
    elsif ($nextPos eq 'U') {
      $seq1 = "-" . $seq1;
      $seq2 = substr($seqs[1],$nextI-1,1) . $seq2;
      $nextI--;
    }
    elsif ($nextPos eq 'L') {
      $seq1 = substr($seqs[0],$nextJ-1,1) . $seq1;
      $seq2 = "-" . $seq2;
      $nextJ--;
    }
    elsif ($nextPos eq 'END') {

      $flagEnd = 1;
    }
    $nextPos = $tb[$nextI][$nextJ];
  }

  return ($seq1, $seq2);
}
